﻿inline void couette(double mason,double posz, double *frep, double *epoten)
{
    ///Fuerzas asociadas al arrastre de Stokes
    *frep=32.0*mason*posz;
    *epoten=16.0*mason*posz*posz;
}

inline void fhs(double dist, double *frep, double *epoten)
{
   // double alpha=3.0e-9;
   double alpha=1.0;
    double c=0.97;
    //double rm6=pow(dist-c,-6.0);
    double rm6=pow(dist,-6.0);
    double rm7=rm6/dist;
    *frep=-alpha*rm6;
    *epoten=alpha/7.0*rm7;
    if(*epoten>100.0)
    {

        cout<<"distancia: "<<dist<<"  "<<*frep<<"  "<<rm6<<"  "<<dist-c<<endl;
        cout<<endl;
    }

}




inline void klingemberg(double dist, double u0, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    if (dist<smed+rcut)
    {
        dummy=u0*exp(-100*(dist-smed));
        ///Para las fuerzas
        *frep=-3.0*dummy;
        ///Para el potencial
        *epoten=dummy;
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void RShiftedLJ(double dist, double sigma, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    double r0,roff,frac2,frac6;
    r0=smed-sigma;
    double rc=smed+2.5*sigma/r0;
    if (dist<rc)
    {

        roff=dist-r0;
        frac2 = (sigma/roff)*(sigma/roff);
        frac6 = frac2*frac2*frac2;

        ///Para las fuerzas
        *frep   = -48.0 * frac6*(frac6 - 0.5) / (roff);
        ///Para el potencial
        *epoten= 4.0*(frac6*frac6-frac6);
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void fbound(particula part, parametros param, double *fpart, double *epoten)
{
    ///Distancia de la particula a las paredes superior e inferior
    ///Pared superior: zpu=+lz/2+a
    ///Pared inferior: zpd=-lz/2-a
	double dzu, dzd;
	double rcut = 0.0;
//	double rcut = param.rcut;
    double alpha=10.0;

    dzu=param.lbox[2]-(part.pos[2]+0.5);
    dzd=part.pos[2]-0.5;

    if(dzu<0.0) {
        ///Interaccion con pared superior
        *fpart=-alpha*(1.-2.*part.pos[2]);
        ///Para el potencial
        *epoten=alpha*(1.-part.pos[2])*part.pos[2];
    }
    if(dzd<0.0){

        ///La particula interactua con la pared inferior
        *fpart=alpha*(1.-2.*part.pos[2]);
        ///Para el potencial
        *epoten=alpha*(1.-part.pos[2])*part.pos[2];
    }
}


inline void magnetic(double dist, double dz, double *fpartxy, double *fpartz, double *epoten)
{
    double Co, Co2;
    double dist3= dist*dist*dist;
    double dist4= dist3*dist;
    double u30dis=1.0/dist4;
    ///Calculamos la fuerza entre iii y jjj
    Co=dz/dist;
	Co2=Co*Co;
	*fpartxy= u30dis*(5.0*Co2-1.0);
	*fpartz= u30dis*(5.0*Co2-3.0);
	*epoten = (1.0-3*Co2)/(dist3);
}

inline double interaccion(particula *part, parametros param)
{

	double dist,dist2,ePoten;
	double fx,fy,fz;
	double dr[3];
	///Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	double fcomp, pcomp;
	///Cutoff a la interaccion
	double rcut=param.rcut;
	///Cutoff interaccion con paredes
	double cutoffwall=param.rcut;
	double eps=1./param.timBrown;
    double alea ;
    double gamma_new;



	///Inicializamos la energia potencial
	ePoten=0.;
	///Inicializamos el tensor de esfuerzos
	double fHS, eHS;
	double fMagxy, fMagz, fMagE;
	fMagxy=fMagz=fMagE=0.0;
	fHS=eHS=0.0;

	///Inicializamos las fuerzas
	for (int iii=0;iii<param.npart;iii++)
	{
	    for(int jjj=0;jjj<3;jjj++)
	    {
	        part[iii].force[jjj]=0.0;
	    }
	}
	///Recorremos todos los pares de particulas
	for(int iii=0;iii<param.npart-1;iii++){
        #ifdef CONFINED
        	///Interaccion de la particula iii con las paredes
            ///Pared superior: zpu=+lz/2+a
            ///Pared inferior: zpd=-lz/2-a
            fbound(part[iii], param, &fcomp, &pcomp);
            part[iii].force[2]+=fcomp;
	    #endif
	   // part[iii].force[0]+=32.0*param.mason*part[iii].pos[2];
	    for(int j=0;j<3;j++)
	    {
            alea=(double(rand())/RAND_MAX-0.5);
            gamma_new=(1.0-eps)*part[iii].gamma+eps*alea;
            part[iii].gamma=gamma_new;
            part[iii].force[j]+=gamma_new*sqrt(2.0*param.timBrown/(3.0*param.lambda))*sqrt((2.0-eps)/eps);
	    }
		for(int jjj=iii+1;jjj<param.npart;jjj++){
		    ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
		    dist = distancia2(part[iii],part[jjj], &param);
		//    if(dist<1.0)
                fhs(dist, &fHS, &eHS);
            // RShiftedLJ(dist, 1.5, 3.0, s1, s2, &fHS, &eHS);
			magnetic(dist,param.dr[2], &fMagxy, &fMagz, &fMagE);

            ///Calculamos la fuerza entre iii y jjj
            fx= (fMagxy+fHS)*(param.dr[0])/dist;
            fy= (fMagxy+fHS)*(param.dr[1])/dist;
            fz= (fMagz+fHS)*(dr[2])/dist;

            #ifdef CALC_STRESS
                sxy+=(param.dr[0])*fy;
                sxz+=(param.dr[0])*fz;
                syz+=(param.dr[1])*fz;
                szz+=(param.dr[2])*fz;
			#endif

			part[iii].force[0]+=fx;
			part[iii].force[1]+=fy;
			part[iii].force[2]+=fz;

			part[jjj].force[0]-=fx;
			part[jjj].force[1]-=fy;
			part[jjj].force[2]-=fz;

			///Acualizamos la energia potencial
			ePoten+=eHS+fMagE;
		}
	}

	#ifdef CONFINED
        ///Consideramos la interaccion con las paredes de la ultima particula
        ///no considerada en el bucle anterior
        fbound(part[param.npart-1], param, &fcomp, &pcomp);
        part[param.npart-1].force[2]+=fcomp;

    #endif
    part[param.npart-1].force[0]+=32.0*param.mason*part[param.npart-1].pos[2];
    for(int j=0;j<3;j++)
    {
        alea=(double(rand())/RAND_MAX-0.5);
        gamma_new=(1.0-eps)*part[param.npart-1].gamma+eps*alea;
        part[param.npart-1].gamma=gamma_new;
        part[param.npart-1].force[j]+=gamma_new*sqrt(2.0*param.timBrown/(3.0*param.lambda))*sqrt((2.0-eps)/eps);
    }
	#ifdef CALC_STRESS
		///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
		spxy=sxy/param.vol;
		spxz=sxz/param.vol;
		spyz=syz/param.vol;
		spzz=szz/param.vol;
	#endif
    //couette(part, param, 1.0, 1.0, param.t);
	return ePoten/param.npart;;
}

inline void idealWallInteraction(particula part, parametros param, double *fpart, double *epoten)
{
    ///Distancia de la particula a las paredes superior e inferior
    ///Pared superior: zpu=+lz/2+a
    ///Pared inferior: zpd=-lz/2-a
	double dzu, dzd;
	double diametro = part.diametro;
	double dummy;

    dzu=param.lbox[2]-(part.pos[2]+0.5*diametro);
    dzd=part.pos[2]-0.5*diametro;

    if(dzu<(0.5*diametro+param.rcut)) {
        ///Interaccion con pared superior
        dummy=param.u0*exp(-100.*(dzu-0.5*diametro));
        ///Para las fuerzas
        *fpart=-3*dummy;
        ///Para el potencial
        *epoten=dummy;
    }
    if(dzd<(diametro+param.rcut)){

        ///La particula interactua con la pared inferior
        dummy=param.u0*exp(-50.*(dzd-0.5*diametro));
        ///Para las fuerzas
        *fpart=3*dummy;
        ///Para el potencial
        *epoten=dummy;
    }
}


inline double interaccionConfined(particula *part, parametros param)
{
    ///Las interacciones consideradas son:
    ///    - Interaccion magnetica entre particulas
    ///    - Interaccion repulsiva particula-particula
    ///    - Interaccion repulsiva particula-pared
	///Si calcStress = true, se calcula la parte correspondiente al potencial
	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
    ///En el calculo del tensor de esfuerzos solamente se tienen en cuenta
    ///la interaccion magnetica y la repulsion particula-particula, pero no
    ///la interaccion particula/caja. COMPROBAR!!!
	///Numero de pares interactuantes
	double dist,dist2,ePoten;
	double fx,fy,fz;
	double dr[3];
	///Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	double fcomp, pcomp;
	//void stress(caja *bx);
	///Cutoff a la interaccion de klingemberg
	double rcutklin=param.rcut;
	///Radio de las particulas interactuantes
	double s1=1.0;
	double s2=1.0;
	///Campo magnetico
	double u0 = param.u0;
	///Cutoff interaccion con paredes
	double cutoffwall=param.rcut;
	///Ponemos las fuerzas de las particulas de la caja a 0
	//bx->initForces();
	///Inicializamos la energia potencial
	ePoten=0.;
	///Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;
	double fHS, eHS;
	double fMagxy, fMagz, fMagE;
	double massi=1.0;
	double massj=1.0;
	double u0Poly;
	fMagxy=fMagz=fMagE=0.0;
	fHS=eHS=0.0;

	///Inicializamos las fuerzas
	for (int iii=0;iii<param.npart;iii++)
	{
	    for(int jjj=0;jjj<3;jjj++)
	    {
	        part[iii].force[jjj]=0.0;
	    }
	}
	///Recorremos todos los pares de particulas
	for(int iii=0;iii<param.npart-1;iii++){
        #ifdef CONFINED
        	///Interaccion de la particula iii con las paredes
            ///Pared superior: zpu=+lz/2+a
            ///Pared inferior: zpd=-lz/2-a
            idealWallInteraction(part[iii], param.u0, param.lbox[2], cutoffwall, &fcomp, &pcomp);
            part[iii].force[2]+=fcomp;
	    #endif
		for(int jjj=iii+1;jjj<param.npart;jjj++){
            s1=part[iii].diametro;
            s2=part[jjj].diametro;
            massi=part[iii].mass;
            massj=part[jjj].mass;
            u0Poly=u0*massi*massj;


		    ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
		    dist = distancia2(part[iii],part[jjj], &param);

            klingemberg(dist, u0Poly, rcutklin, s1, s2, &fHS, &eHS);
            // RShiftedLJ(dist, 1.5, 3.0, s1, s2, &fHS, &eHS);
			//magnetic(dist,param.dr[2], &fMagxy, &fMagz, &fMagE);


            ///Calculamos la fuerza entre iii y jjj
            fx= (fMagxy+fHS)*(param.dr[0])/dist;
            fy= (fMagxy+fHS)*(param.dr[1])/dist;
            fz= (fMagz+fHS)*(dr[2])/dist;

//            if(fx>10. || fy>10. || fz>10.)
//            {
//                cout<<fHS<<"  "<<dist<<endl;
//                cout<<endl;
//            }

            #ifdef CALC_STRESS
                sxy+=(param.dr[0])*fy;
                sxz+=(param.dr[0])*fz;
                syz+=(param.dr[1])*fz;
                szz+=(param.dr[2])*fz;
			#endif

			part[iii].force[0]+=fx;
			part[iii].force[1]+=fy;
			part[iii].force[2]+=fz;

			part[jjj].force[0]-=fx;
			part[jjj].force[1]-=fy;
			part[jjj].force[2]-=fz;

			///Acualizamos la energia potencial
			ePoten+=eHS+fMagE;
		}
	}
	//cout<<"***************+  "<<ePoten<<endl;
	#ifdef CONFINED
        ///Consideramos la interaccion con las paredes de la ultima particula
        ///no considerada en el bucle anterior
         idealWallInteraction(part[bx->npart-1], param, &fcomp, &pcomp);
        part[npart-1].force[2]+=fcomp;
    #endif
	#ifdef CALC_STRESS
		///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
		spxy=sxy/param.vol;
		spxz=sxz/param.vol;
		spyz=syz/param.vol;
		spzz=szz/param.vol;
	#endif
//    couette(part, param, 1.0, 1.0, param.t);
	return ePoten/param.npart;;
}

//inline double interaccionConfined(caja *bx)
//{
//    ///Las interacciones consideradas son:
//    ///    - Interaccion magnetica entre particulas
//    ///    - Interaccion repulsiva particula-particula
//    ///    - Interaccion repulsiva particula-pared
//	///Si calcStress = true, se calcula la parte correspondiente al potencial
//	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
//    ///En el calculo del tensor de esfuerzos solamente se tienen en cuenta
//    ///la interaccion magnetica y la repulsion particula-particula, pero no
//    ///la interaccion particula/caja. COMPROBAR!!!
//	///Numero de pares interactuantes
//	double dist,dist2,ePoten;
//	double fx,fy,fz;
//	double virial=0;
//	///Tensor de esfuerzos
//	double sxy, sxz, syz, szz;
//	double fcomp, pcomp;
//	void stress(caja *bx);
//	///Cutoff a la interaccion de klingemberg
//	double rcutklin=bx->param.rcut;
//	///Radio de las particulas interactuantes
//	double s1=1.0;
//	double s2=1.0;
//	///Campo magnetico
//	double u0 = bx->u0;
//	///Cutoff interaccion con paredes
//	double cutoffwall=bx->param.rcut;
//	///Ponemos las fuerzas de las particulas de la caja a 0
//	//bx->initForces();
//	///Inicializamos la energia potencial
//	ePoten=0.;
//	///Inicializamos el tensor de esfuerzos
//	sxy=sxz=syz=szz=0.0;
//	double fHS, eHS;
//	double fMagxy, fMagz, fMagE;
//	double massi=1.0;
//	double massj=1.0;
//	double u0Poly;
//	double dr[3];
//	int np=bx->npart;
//	double lx,ly,lz;
//	lx=bx->lbox[0];
//	ly=bx->lbox[1];
//	lz=bx->lbox[2];
//	///Recorremos todos los pares de particulas
//	for(int iii=0;iii<np-1;iii++){
//	    ///Interaccion de la particula iii con las paredes
//	    ///Pared superior: zpu=+lz/2+a
//	    ///Pared inferior: zpd=-lz/2-a
//
//	    idealWallInteraction(bx->part[iii], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);
//	    bx->part[iii].force[2]+=fcomp;
//		for(int jjj=iii+1;jjj<np;jjj++){
//            #ifdef POLYDISPERSE
//			s1=bx->part[iii].radio;
//			s2=bx->part[jjj].radio;
//			massi=bx->part[iii].mass;
//			massj=bx->part[jjj].mass;
//			u0Poly=u0*massi*massj;
//			#endif
////
////		    ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
////    dr[0]=bx->part[jjj].pos[0]-bx->part[iii].pos[0];
////    dr[0]-=lx*anint(dr[0]/lx); //Entero + proximo manteniendo signo
////    /*Y*/
////    dr[1]=bx->part[jjj].pos[1]-bx->part[iii].pos[1];
////    dr[1]-=ly*anint(dr[1]/ly); //Entero + proximo manteniendo signo
////    /*Z*/
////    dr[2]=bx->part[jjj].pos[2]-bx->part[iii].pos[2];
////    #ifndef CONFINED
////    dr[2]-=lz*anint(dr[2]/lz); //Entero + proximo manteniendo signo
////    #endif
////
////	dist2=dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2];
//
//			dist2= bx->distancia2(iii,jjj);
//
//			//dist2=4.;
//			dist=sqrt(dist2);
//
//			#ifdef POLYDISPERSE
//            klingemberg(dist, u0Poly, rcutklin, s1, s2, &fHS, &eHS);
//			magnetic(dist,bx->dr[2], u0Poly, &fMagxy, &fMagz, &fMagE);
//			//fMagxy=fMagz=fMagE=0.0;
//			//fMagxy=fMagz=fMagE=fHS=eHS=0.0;
//            #endif
//
//            #ifndef POLYDISPERSE
//			klingemberg(dist, u0, rcutklin, s1, s2, &fHS, &eHS);
//			magnetic(dist,bx->dr[2], u0, &fMagxy, &fMagz, &fMagE);
//			//fMagxy=fMagz=fMagE=0.0;
//			#endif
//
//            ///Calculamos la fuerza entre iii y jjj
//            fx= (fMagxy+fHS)*(bx->dr[0])/dist;
//            fy= (fMagxy+fHS)*(bx->dr[1])/dist;
//            fz= (fMagz+fHS)*(bx->dr[2])/dist;
//
//            #ifdef CALC_STRESS
//                sxy+=(dr[0])*fy;
//                sxz+=(dr[0])*fz;
//                syz+=(dr[1])*fz;
//                szz+=(dr[2])*fz;
//			#endif
//
//			bx->part[iii].force[0]+=fx;
//			bx->part[iii].force[1]+=fy;
//			bx->part[iii].force[2]+=fz;
//
//			bx->part[jjj].force[0]-=fx;
//			bx->part[jjj].force[1]-=fy;
//			bx->part[jjj].force[2]-=fz;
//
//			///Acualizamos la energia potencial
//			ePoten+=eHS+fMagE;
////			///Actualizamos el virial
////			virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);
//
//		}
//	}
//	///Consideramos la interaccion con las paredes de la ultima particula
//	///no considerada en el bucle anterior
//	 idealWallInteraction(bx->part[bx->npart-1], u0, bx->lbox[2], cutoffwall, &fcomp, &pcomp);
//
//    bx->part[np-1].force[2]+=fcomp;
//;
//	#ifdef CALC_STRESS
//		///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
//		bx->spxy=sxy/bx->vol;
//		bx->spxz=sxz/bx->vol;
//		bx->spyz=syz/bx->vol;
//		bx->spzz=szz/bx->vol;
//
//	#endif
//
//	bx->epoten=ePoten/np;
//	bx->virial=virial;
//	return ePoten;
//}

